!***********************************************************
!
program plane_2_line
!BE CAREFUL WITH FORTRAN COMPILER
!USE OPTION ifort -assume byterecl
!***********************************************************
USE ir_precision
USE lib_vtk_io

implicit none

 ! Variables reading from Incompact3D, real(8) for double precision
 integer,parameter :: nx=1,ny=65,nz=64, new_nz=1
 real(8), parameter :: lx=31.4, ly=2., lz=6.28
 real(8),dimension(nx,ny,nz) :: read_data
 real(8),dimension(ny) :: yp, write_data1
 real(8) :: dx, dy, dz
 integer :: i,j,k,count
 
 ! Variables writing to Paraview
 integer, parameter :: step=1
 integer(I4P):: e_io
 integer(I4P), parameter :: nx1=int(1,i4p),nx2=int(nx/step,i4p)
 integer(I4P), parameter :: ny1=int(1,i4p),ny2=int(1+(ny-1)/step,i4p)
 integer(I4P), parameter :: nz1=int(1,i4p),nz2=int(new_nz/step,i4p)
 integer(I4P), parameter :: nn=(nx2-nx1+1)*(ny2-ny1+1)*(nz2-nz1+1)
 real(R4P), dimension(nx1:nx2,ny1:ny2,nz1:nz2) :: x,y,z
 real(R4P), dimension(nx1:nx2,ny1:ny2,nz1:nz2) :: write_data

 ! Variables pour gestion liste fichiers
 integer :: nca, l
 character(20):: nf, line_file

 ! Récupération nombre arguments (Variables for managing files list)
 nf = ''
 nca = command_argument_count()
 if (nca==0) then
  write(*,*) 'Programme quits. No files to read / write ...'
  stop
 else
  write(*,*) nca,' files to read and export for Paraview'
 endif
 
 ! dx et dz fonction des CL (dx and dz function of CL)
 dx=lx/(nx)!nx ! nx pour périodique, nx-1 sinon (for periodic dx=lx/(nx-1), for non-periodic dx=lx/(nx) )
 dz=lz/(nz)!nz ! nz pour périodique, nz-1 sinon (for periodic dz=dz/(nz-1), for non-periodic dz=lz/(nz) )
 
 ! Reading the mesh
 write(*,*) 'Reading the grids y'
 open(21,file='yp.dat',form='formatted')
 do j=1,ny
  read(21,*) i,yp(j)
 enddo
 close(21)
 write(*,*) 'Grids are read'
 
 ! Calculating coordinates for Paraview
 DO k=nz1,nz2
  DO j=ny1,ny2
   DO i=nx1,nx2
    x(i,j,k)=real(   ((i-1)*step)*dx, R4P)
    y(i,j,k)=real( yp(1+(j-1)*step), R4P)
    z(i,j,k)=real(   ((k-1)*step)*dz, R4P)
   ENDDO
  ENDDO
 ENDDO

 
 ! Database file for structured grid
 write(*,*) 'Initializing file for Paraview'
 e_io = vtk_ini_xml(output_format='raw',filename='XML_STRG.vts',mesh_topology='StructuredGrid',&
 nx1=nx1,nx2=nx2,ny1=ny1,ny2=ny2,nz1=nz1,nz2=nz2)
 ! Geometry file
 write(*,*) 'Initializing geometry'
 e_io = vtk_geo_xml(nx1=nx1,nx2=nx2,ny1=ny1,ny2=ny2,nz1=nz1,nz2=nz2,nn=nn,x=x,y=y,z=z)
 ! Opening data file for Paraview
 write(*,*) 'Initializing data'
 e_io = vtk_dat_xml(var_location='node',var_block_action='open')

 ! Loop playback file + Interpolation + Writing for paraview
 ! ATTENTION: SI DATA EN PROVENANCE IDRIS ALORS CONVERT='big_endian' (IF DATA FROM IDRIS THEN CONVERT)

 do l=1,nca
  call get_command_argument(l,nf)
  ! Reading file
  write(*,*) 'Reading file ',trim(nf)
  OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
     ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  COUNT = 1
  DO k=1,nz
   DO j=1,ny
    DO i=1,nx
     READ(21,REC=COUNT) read_data(i,j,k)
     COUNT = COUNT + 1
    ENDDO
   ENDDO
  ENDDO
  CLOSE(21)
  write(*,*) 'Reading file ',trim(nf),' completed'
  write(*,*) 'Max & Min : ',maxval(read_data),minval(read_data)

  write_data1=0.
  !Plane to Line
  DO k=1,nz
   DO j=1,ny
    DO i=1,nx
     write_data1(j) = write_data1(j)+read_data(i,j,k)/nz
    ENDDO
   ENDDO
  ENDDO
  
  line_file=trim(nf)//'_line'

  ! Writing file
  write(*,*) 'Writing file ',line_file
  OPEN(21,FILE=line_file,FORM='UNFORMATTED',&
     ACCESS='DIRECT', RECL=8)!, CONVERT='big_endian')
  COUNT = 1
  DO j=1,ny
    write(21,REC=j) write_data1(j)
    COUNT = COUNT + 1
  ENDDO
  CLOSE(21)
  
  ! Interpolation
  write(*,*) 'Start interpolation'
  DO k=nz1,nz2
   DO j=ny1,ny2
    DO i=nx1,nx2
     write_data(i,j,k) = real( write_data1(1+(j-1)*step), R4P)
    ENDDO
   ENDDO
  ENDDO
  write(*,*) 'Interpolating completed'
  ! Writing Paraview
  write(*,*) 'Adding file ',line_file,' in Paraview'
  e_io = vtk_var_xml(nc_nn=nn,varname=line_file,var=write_data)
  write(*,*) 'File ',line_file,' added to Paraview'
 ENDDO
 ! end loop
 
 ! Close all Paraview files
 write(*,*) 'Closing Paraview files (*.vts)'
 e_io = vtk_dat_xml(var_location='node',var_block_action='close')
 e_io = vtk_geo_xml()
 e_io = vtk_end_xml()


 write(*,*) 'Terminus'

end program plane_2_line

